%% Set up everything that doesn't need to be done in the loop
% Load Excel data with ethanol plant capacities and locations
ethanolPlants = readtable("../Ethanol Plant Data/Plant Effective Radius.xlsx");

% Extract plant locations as row vectors
plantLatLong = [ethanolPlants.lat ethanolPlants.lng]';

% Years of analysis
T = 2001:2018;
% Variable names for the ethanol index
varNames1 = "Z1_"+T;
varNames2 = "Z1CF_"+T;
varNames3 = "Z2_"+T;
varNames4 = "Z2CF_"+T;

% Radius of the earth in miles; needed for distance calc
R = 3959;

%% Load field centroid data
% Load Excel data with sample pixels
pixels = readtable("../CSB/fieldCentroids1.csv");

%% Analyze data
% Number of pixels in pixel data
n_pix = height(pixels);

pixLatLong = [pixels.lat pixels.lng];

%% Compute (n_pix x n_plants) distance matrix
% Cell (i,j) is distance between pixel i and plant j

% Define constants for the haversine calculation
phi1 = pixLatLong(:,1)*pi()/180; % Pixel latitude in radians
phi2 = plantLatLong(1,:)*pi()/180; % Plant latitude in radians
deltaPhi = phi1-phi2; % Difference between Pixel lat and Plant lat in radians
deltaLambda = (pixLatLong(:,2)-plantLatLong(2,:))*pi()/180; % Difference between Pixel long and Plant long in radians

a = sin(deltaPhi/2).*sin(deltaPhi/2) + cos(phi1).*cos(phi2).*sin(deltaLambda/2).*sin(deltaLambda/2);
c = 2*atan2(sqrt(a),sqrt(1-a));

distances = R*c; % Distance measured in miles

%% Loop over years and calculate the distance index for each pixel

% Create place holder for the distance index
% One value for each year for each pixel
index1 = zeros(n_pix,length(T));
indexCF1 = index1;
index2 = index1;
indexCF2 = index1;

for y=1:length(T)
    % Compute ethanol index with index(y,index_string,[params]) function
    index1(:,y) = index(y,"d1_",T,ethanolPlants,distances);
    indexCF1(:,y) = index(y,"dCF1_",T,ethanolPlants,distances);
    index2(:,y) = index(y,"d2_",T,ethanolPlants,distances);
    indexCF2(:,y) = index(y,"dCF2_",T,ethanolPlants,distances);
end

%% Export data
% Convert distanceIndex matrix to table
T1 = array2table(index1,"VariableNames",varNames1);
T2 = array2table(indexCF1,"VariableNames",varNames2);
T3 = array2table(index2,"VariableNames",varNames3);
T4 = array2table(indexCF2,"VariableNames",varNames4);

exportTable = [pixels T1 T2 T3 T4];
writetable(exportTable,"../Excel files/Matlab Output/datafile.csv");
    
%% Function - Compute ethanol demand index
function Z = index(y,index_string,T,ethanolPlants,distances)
    % Get current year
    yr = T(y);
    
    % Use year to get the variable name of the d_bar column we are looking for
    var_string = index_string+num2str(yr);

    % Get the d_bar column as a row vector
    d_bar = (ethanolPlants.(var_string))';

    % Compute index = max value of d_bar_j - d_ij across plants j, or 0 if
    % negative
    Z = max(max(d_bar-distances,[],2),0);
end
